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CHAPTER 1 

SUPPLY AND PRODUCTION NETWORKS: FROM THE 
BULLWHIP EFFECT TO BUSINESS CYCLES 
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Institute for Transport & Economics, Dresden University of Technology 
Andreas- Schubert-Str. 23, D-01062 Dresden, Germany 
E-mail: helbing@trafficforum.org, trajfic@stefanlaemmer.de 



Network theory is rapidly changing our understanding of complex sys- 
tems, but the relevance of topological features for the dynamic behavior 
of metabolic networks, food webs, production systems, information net- 
works, or cascade failures of power grids remains to be explored. Based 
on a simple model of supply networks, we offer an interpretation of in- 
stabilities and oscillations observed in biological, ecological, economic, 
and engineering systems. 

We find that most supply networks display damped oscillations, even 
when their units - and linear chains of these units - behave in a non- 
oscillatory way. Moreover, networks of damped oscillators tend to pro- 
duce growing oscillations. This surprising behavior offers, for example, 
a new interpretation of business cycles and of oscillating or pulsating 
processes. The network structure of material flows itself turns out to be 
a source of instability, and cyclical variations are an inherent feature of 
decentralized adjustments. In particular, we show how to treat produc- 
tion and supply networks as transport problems governed by balance 
equations and equations for the adaptation of production speeds. 

The stability and dynamic behavior of supply networks is investi- 
gated for different topologies, including sequential supply chains, "supply 
circles", "supply ladders", and "supply hierarchies". Moreover, analyti- 
cal conditions for absolute and convective instabilities are derived. The 
empirically observed bullwhip effect in supply chains is explained as a 
form of convective instability based on resonance effects. An applica- 
tion of this theory to the optimization of production networks has large 
optimization potentials. 
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1. Introduction 



Supply chain management is a major subject in economics, as it signifi- 
cantly influences the efficiency of production processes. Many related stud- 
ies focus on subjects such as optimum buffer sizes and stock levels! I I I I 
However, the optimalstructure of supply and distribution networks is also 
an important topicPl w hich connects this s cientific field with the statis- 
tical physics of networks-EEMinillllllini p 

roblems like the adaptivity 
and robustness of supply networks have not yet been well covered and are 
certainly an important field for future research. 

In this chapter, we will focus on the dynamical properties and lin- 
ear stability of supply networks in dependence of the network topology 
or, in other words, the supply matrix (input matrix). Presently, there are 
only a few results on this subject, since the response of supply networks 
to varying demand is not a trivial problem, as we will see. Some "fl uid 
models" have, however, been proposed to study this subject. Daganzo, - ^ 
for example, applies the method of cumulative counts—^ and a variant 
of the Lighthill-Whitham traffic model to study the so-called bullwhip 
effect ! 1 7 l l« I ^ I ^Pl | 22 | 2:j | 24 | 25|2tj | 27 EHlThis effect has, for example, been 
reported for beer distribution^ - iHI k u t similar dynamical effects are also 
known for other distribution or transportation chains. It describes increas- 
ing oscillations in the delivery rates and stock levels from one supplier to 
the next upstream supplier delivering to him. _______ 

Similar models are studied by Armbruster et a/ P —————31 Most pub- 
lications, however, do not investigate the impact of the netw ork to pology, 
but focus on sequential supply chains or re-entrant problems An ex- 
ception is the transfe r fu nction approach by Dejonckheere et o/!— — ' as well 
as Disney and Towill) - 1 Ponzi, Yasutomi, and KanekJ - 1 have coupled a 
supply chain model with a model of price dynamics, while Witt and Surl— ^ 
have suggested to model business cycles analogously to stop-and-gotraf- 
fic. Similar suggestions have been made by Daganzo>23 and HelbingJjpWe 
should also note the relationship with driven many-particle models E_ This 
is the basis of event-driven simulations of production systems, which are 
often used in logistics software. 

In a previous publication, one of the authors has proposed a rather gen- 
eral model for suppl y n etworks , and has connected it to queueing theory 
and macroeconomics!— 21 That paper also presents numerical studies of the 
bullwhip effect in sequential supply chains and the effect of heterogeneous 
production units on the dynamics of supply networks. Moreover, the de- 
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pcndcncc of the maximum oscillation amplitude on the model parameters 
has been numerically studied, in particular the phase transition from stable 
to bullwhip behavior. The underlying simulation models are of non-linear 
naturep^l so that p hen omena such as wave selection and synchronization 
have been observed!^ Finally, a large variety of management strategies 
includin g the e ffect of forcasting inventories has been studied together with 
NagatanP^E] The latter has also investigated the dynamical effect of mul- 
tiple production lines@3 

In the following review, we will mainly focus on a linearized model of 
supply networks, which could be view ed as a dynamical versio 

rP 

of Leon- 

tief's classical input-output modelE^It is expected to be valid close to the 
stationary state of the supply system or, in other words, close to the (com- 
monly assumed) economic equilibrium. This model allows to study many 
properties of supply networks in an analytical way, in particular the effects 
of delayed adaptation of production rates or of a forecasting of inventories. 
In this way, we will be able to understand many numerically and empiri- 
cally observed effects, such as resonance effects and the stability properties 
of supply networks in dependence of the network topology. 

Our review is structured as follows: Section [2] introduces our model of 
supply networks and linearizes it around the stationary state. Section 12.21 
discusses the bullwhip effect for sequential supply chains. Surprisingly, the 
amplification of amplitudes does not require unstable eigenvalues, but it is 
based on a resonance effect. Sec. lH.II presents general methods of solution for 
arbitrary topologies of supply networks. It reveals some useful properties of 
the eigenvalues characterizing the stability of supply systems. Interestingly 
enough, it turns out that all supply networks can be mapped on formulas for 
linear (serial) supply chains. Sec. 13.21 will investigate, how the dynamic be- 
havior of supply networks depends on their network structure, while Sec.0] 
gives an interpretation of business cycles based on material flow networks. 
Sec. El finally summarizes our results and points to some future research 
directions. 

2. Input-Output Model of Supply Networks 

Our production model assumes u production units or suppliers j £ 
{1, . . . , u} which deliver d^ products of kind i 6 {l,...,p} per production 
cycle to other suppliers and consume Ckj goods of kind k per production 
cycle. The coefficients Ckj and are determined by the respective produc- 
tion process, and the number of production cycles per unit time (e.g. per 
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4 Dirk Helbing and Stefan hammer 

day) is given by the production speed Qj{t). That is, supplier j requires an 
average time interval of 1/Qj(t) to produce and deliver dij units of good 
i. The temporal change in the number iVj(t) of goods of kind i available in 
the system is given by the difference between the inflow 

u 
3=0 

and the outflow 

u+l 

Qr t (*) = E c ^^(*)- ( 2 ) 

In other words, it is determined by the overall production rates dijQj(t) of 
all suppliers j minus their overall consumption rates CijQj(t): 



"' V Qt{t)-QT\t) = Y,diiQj{t) 



dt 

3 = 1 



3 = 1 



(3) 



supply 



In this dynamic variant of Leontief 's classical input-output modelp^ the 
quantity 

Yi(t)= Cj.u+iQn+ift) - cf«QoW (4) 

consumption and losses inflow of resources 

comprises the consumption rate of goods i, losses, and waste (the "export" 
of material), minus the inflows into the considered system (the "imports"). 
In the following, we will assume that the quantities are measured in a way 
that < Cij , dij < 1 (for 1 < i < p, 1 < j < u) and that the "normalization 
conditions" 



dio = 1 - y ]djj > , c i>u+ i = 1 - 2J <Hj > (5) 

3=1 3=1 

are fulfilled. Equations (0 can then be interpreted as conservation equations 
for the flows of goods. 

2.1. Adaptation of Production Speeds 

In addition, we have formulated an equation for the production or delivery 
rate Qjit). Changes in the consumption rate Yi(t) sooner or later require 
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Supply Networks: From the Bullwhip Effect to Business Cycles 5 

an adaptation of Qj(t). For the temporal change dQj/dt of the delivery 
rate we will assume: 

^l=F j ({N i (t)}ddN i /dt},{Q l (t)}). (6) 

Herein, the curly brackets indicate that the so-called management or control 
function Fj{. . . ) may depend on all inventories Ni(t) with i £ {1, . . . , p}, 
their derivatives dNi/dt, and/or all production speeds Qi(t) with I £ 
{1, . . . , u}. Some reasonable specifications will be discussed later on. 

2.2. Modelling Sequential Supply Chains 



Fig. 1. Illustration of the linear supply chain treated in this chapter, including the key 
variables of the model. Circles represent different suppliers i, JVj their respective stock 
levels, and Qi the delivery rate to supplier i or the production speed of this supplier. 
i = corresponds to the resource sector generating the basic products and i = u + 1 to 
the consumers. 



For simplicity, let us first investigate a model of sequential supply chains 
(see Fig.^), which corresponds to djj = 5™ and = S^ij. In other words, 
the products i are directly associated with producers j and we have an input 
matrix of the form 

fO 1 0\ 
10 

1 0. (7) 
1 
Vo 0/ 

This implies Qf{t) = Q t (t) = Q^\(t) and Q° ut (t) = Q l+1 (t) = Qf +1 (t), so 
that the inventory of goods of kind i changes in time t according to 

= Q?(t) - or (*) = Qi® - Qi+i(t) ■ (8) 

The assumed model consists of a series of u suppliers i, which receive prod- 
ucts of kind i — 1 from the next "upstream" supplier i — 1 and generate prod- 
ucts of kind i for the next "downstream" supplier i + 1 at a rate Qi(t) 
The final products are delivered at the rate Q u (t) and removed from the 
system with the consumption rate Y u (t) = Q u +\(t). 

In the following, we give some possible specifications of equation © for 
the production rates: 
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The delivery rate Qi(t) may, for example, be adapted to a certain 
desired rate Wi (Ni, dNi/dt) according to the equation 



dQj 
~df 



1 



m 

dt 



Qi(t) 



(9) 



where denotes the relaxation time. 
A special case of this is the equation 



dQj 
dt 



1 
1 



(Ni(t)) 



w, 



Wi[Ni(t + At) 



dWi(Ni) M \dN,, 



dN 



dt 



-Qi(t) 



=/3» 
i(t) 



(10) 



Herein, the desired production speed Wi(Ni) is monotonously 
falling with increasing inventories iV», which are forecasted over 
a time period At. 

A further specification is given by the equation 



dQj 
~~dT 



1 



N? ~ Ni(t) 



(11) 



which assumes that the production rate is controlled in a way that 
tries to reach some optimal inventory N® and production rate Q®, 
and attempts to miminize changes dN / dt in the inventory to reach 
a constant work in progress (CONWIP strategy). Note that at least 
one of the parameters Tj, Tj, /3j and could be dropped. 
• In other cases, it can be reasonable to work with a model focussing 
on relative changes in the variables: 



dQi/dt 
ft 



N(t) 



1 



dNi/dt 
N(t) 



QL 



l 



(12) 



This model assumes large production rates when the inventory is 
low and prevents that Qi(t) can fall below zero. Apart from this, 
Ni(t) and Qi(t) are adjusted to some values N® and which 
are desireable from a production perspective (i.e. in order to cope 
with the breakdown of machines, variations in the consumption 
rate, etc.) Moreover, the control strategy (|12[) counteracts temporal 
changes dNi j dt in the inventory. For comparison with the previous 
strategies, we set Tj = 1/Q", n = N®/Di, — fii/Nf, and = 
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Supply Networks: From the Bullwhip Effect to Business Cycles 7 

The above control strategies appear to be appropriate to keep the inven- 
tories Ni(t) stationary, to maintain a certain optimal inventory JV] 5 (in 
order to cope with stochastic variations due to machine breakdowns etc.), 
and to operate with the equilibrium production rates Q° = Y°, where Y® 
denotes the average consumption rate. However, the consumption rate is 
typically subject to perturbations, which may cause a bullwhip effect, i.e. 
growing variations in the stock levels and deliveries of upstream suppliers 
(see Sec. ED . 

2.3. More Detailed Derivation of the Production Dynamics 
(with Dieter Armbruster) 

Let us assume a sequential supply chain and let Qi(t) be the influx in the 
last time period T = 24 hours from the supplier i — 1 to producer i. This 
influx is equal to the release J2j_i(t) over the last 24 hours: 

Qi(t) = Ri-i(t). (13) 

Based on the current state (the inventory and fluxes at time t and before), 
a price Pi{t) is set according to some pricing policy and instantaneously 
communicated downstream to the customer i + 1. Then, based on price, 
fluxes, and inventory, the order quantity Di + \{t) is decided at the end of 
day t according to some order policy Wi+i: 

D i+1 {t)=W i+l (N i+ x{t),P i {t),...') . (14) 

It determines the upstream release Ri(t + T) on the next day t + T = t + 1: 

Ri(t + T)=Ri(t + l)=D i+1 (t). (15) 
Altogether, this implies 

Qi(t + T) = Qi(t + 1) = ifc_i(t + 1) = Di(t) = Wi (Ni(t), . . . ) (16) 
and the following equation for the change in the production rate: 

d& w Aftw = i m+T) _ QM = i [ Wi ( m)t . ..y Qiit )\ . a?) 

This is consistent with Eq. ©. Moreover, we obtain the usual balance 
equation for the change of inventory in time: 

dN, AN t Ni(t + T)-Ni(t) 

-27--^ = — f — = N ^ + !) - N *® = - 

= Qi(t) - Q i+ i(t) . (18) 



February 2, 2008 11:33 



WSPC/Trim Size: 9in x 6in for Review Volume 



hclbinglaemmer 



Dirk Helbing and Stefan hammer 



2.4. Dynamic Solution and Resonance Effects (with 
Tadeusz Platkowski) 

Let us now calculate the dynamic solution of the sequential supply chain 
model for cases where the valu es iV and Q° correspond to the stationary 
state of the production system.lSESl Then, the linearized model equations 
for the control approaches JSJ to (|12fl are exactly the same. Representing 
the deviation of the inventory from the stationary one by Ui(t) = Ni(t) — N® 
and the deviation of the delivery rate by qi(t) = Qi(t) — Q®, they read 



dt 



1 

T, 



rn(t) dn t ' 

Pt-TT- Wi\t) 

Ti dt 



(19) 



Moreover, the linearized equations for the inventories are given by 

— = qi{t) - q l+ i(t) . 

Deriving Eq. (|19|l with respect to t and inserting Eq. (|2U[1 results in the 
following set of second-order differential equations: 

d 2 q l [fii + ei) dqi 1 1 
dt 2 




Ti 



Qi+iit) | dq i+ i 



dt 



(20) 



(21) 



=/<(*) 



This corresponds to the differential equation of a damped harmonic oscil- 
lator with damping constant 7, eigenfrequency u>, and driving term fi(t). 
The eigenvalues of these equations are 



Ai,2 = -7< ± V 



.1 — 



1 

2T, 



(Pi + £i) T v 7 ^ + " 



(22) 



For (Pi + e,) > their real parts are always negative, corresponding to a 
stable behavior in time. Nevertheless, we will find a convective instability, 
i.e. the oscillation amplitude can grow from one supplier to the next one 
upstream. 

Assuming periodic oscillations of the form fi(t) = /? cos(crf), the general 
solution of Eq. (|21|l is of the form 

q t (t) = f^F, cos(at + + D?e-^* cos(ad + 9,) , (23) 

where the parameters and 8i depend on the initial conditions. The other 
parameters are given by 

27;a a(Pi + e t ) 



a 2 T t - l/ n 



(24) 
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n* = vV - 7« 2 = 7f V W - (ft + ei ) 2 / 4 , (25) 



and 



1 Ts 

Fl = ^/W^W^W^ ~ ^Wt^TJWT^w+W ' (26) 

The dependence on the eigenfrequency u>i is important for understanding 
the occuring resonance effect, which is particularly likely to appear, if the 
oscillation frequency a of the consumption rate is close to one of the reso- 
nance frequencies u>i. After a transient time much longer than 1/7,; we find 



qi(t) = f°FiCOs(at + ipi). (27) 
Equations l(2Tjl and (|T7J> imply 

= /°_ 1 cos(ai + ^ + <S l ) (28) 



r« etc 



with 



tan^ = aftri and = fi^VOTn? + ("A) 2 • (29) 

Therefore, the set of equations (|21|l can be solved successively, starting with 
i = u and progressing to lower values of i. 



2.5. T/ie Bullwhip Effect 



Fig. 2. Frequency response for different ft and Tj = Tj = ej = 1. For small ft, corre- 
sponding to a small prognosis time horizon At, a resonance effect with an amplification 
factor greater than 1 can be observed. Perturbations with a frequency a close to the 
cigcnfrcquencies uii = l/yTV^ are amplified and cause variations in stock levels and 
deliveries to grow along the supply chain. This is responsible for the bullwhip effect. 



The oscillation amplitude increases from one supplier to the next up- 
stream one, if 

fi-i _ f x + a 2 [er{e l + 2p i )-2T l /T l ]+a i T? \- 1/2 > % _ 



f? I (lM) 2 + (aft) 2 
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One can see that this resonance effect can occur for < a 2 < 2/(TjTi) — 
ej(e, + 2(3i) /Tf . Therefore, variations in the consumption rate are magnified 
under the instability condition 

T i >e i T i {/3 i + e i /2). (31) 

Supply chains show the bullwhip effect (which corresponds to the phe- 
nomenon of convective, i.e. upstream amplification), if the adaptation time 
Ti is too large, if there is no adaptation to some equilibrium production 
speed (corresponding to ej = 0), or if the production management reacts 
too strong to deviations of the actual stock level Ni from the desired one 
JV? (corresponding to a small value of Ti), see Fig. The latter is very 
surprising, as it implies that the strategy 

^T = 7j!-[N?-m)}, (32) 
at Tin 

which tries to maintain a constant work in progress iVj(t) = N®, would 
ultimately lead to an undesireable bullwhip effect, given that production 
units are adjusted individually, i.e. in a decentralized way. In contrast, the 
management strategy 

f 4{- ftf f -*(«>l} < 33 > 

would avoid this problem, but it would not maintain a constant work in 
progress. 



Fig. 3. (a) Plot of the maximum amplitude of oscillation in the inventories as a function 
of the adaptation time T. (b) Plot of the maximum amplitude of the inve ntor ies as a 
function of the time horizon At for an adaptation time of T = 2. (After RefE3| 



The control strategy (|32|) with a sufficiently large value of Tj would fulfill 
both requirements. Having the approximate relation (|10|l in mind, a forecast 
with a sufficiently long prognosis time horizon At (implying large values 
of Pi) is favourable for production stability. Delays in the determination 
of the inventories Ni, corresponding to negative values of At and /3j, are 
destabilizing. Note that the values of At which are sufficient to stabilize the 
system are often much smaller than the adaptation time Tj.ED 
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3. Network Effects (with Peter Seba) 

The question is now, whether and how the bullwhip effect can be gener- 
alized to supply networks and how the dynamic behavior depends on the 
respective network structure. Linearization of the adaptation equation © 
leads to 

= - E *w«) - E w ^ - E • ( 34 ) 

ill 

However, in order to avoid a vast number of parameters Vji, Wji, Xji and 
to gain analytical results, we will assume product- and sector-independent 
parameters Vjk — VSjk, Wjk — WSjk, and Xji = Sji in the following. This 
case corresponds to the situation that each production unit j is character- 
ized by one dominating product i = j, which also dominates production 
control. For this reason, we additionally set dy = Sij. Generalizations are 
discussed later (see Sec. l7.4f) . 

Using vector notation, we can write the resulting system of differential 
equations as 

^=Sq(t)-v(t) (35) 



and 



^ = -VH(t)-W^-q{t)- (36) 



Inserting Eq. l|3l))) into Eq. |JSSJ gives 

dt \ q J \q J \ Wy 

with 



= M ( ■;) + ( J*) (37) 



m -Ue!-e-Vs) ^ 

and the supply matrix S = D — C = (dij — Cij ) . For comparison of Eq. (|19(l 
the one above, one has to scale the time by introducing the unit time T/e 
and to set V = l/(re), W = 0/e. 

3.1. General Methods of Solution 

It is possible to rewrite the system of 2u first order differential equations 
(137(1 in the form of a system of u differential equations of second order: 

|| + (E + WS)^ + VSq\t) = g(t) , (39) 
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where 



m = vm + wf t . 



Introducing the Fourier transforms 
F(a) - ' 



and 



G(a) 



'2tt 



dt e~ iat q(t) 



dt e- iat g(t) 



reduces the problem to solving 

[-a 2 + ia(E + WS) + VS] F(a) 



G{a) 



(40) 



(41) 



(42) 



(43) 



with given G{a). The variable a represents the perturbation frequencies, 
and the general solution of Eq. (|39l) is given by 



da e lat F(a) 



(44) 



Note that there exists a matrix T which allows one to transform the 
matrix E — S via T _1 (E — S)T = J into either a diagonal or a Jordan 
normal form J. Defining x(r) = T _1 <f(r) and h(t) = T _1 g(t), we obtain 
the coupled set of second-order differential equations 



— + 2^—+^ Xl (t) 



Wx l+ i{t) + V 



dx 



i+i 



dt 



+ lH(t), (45) 



where 

7i = [1 + W(l - J«)]/2 , Ui = [V(l - J«)] 1/2 and ./,, . , . (46) 

This can be interpreted as a set of equations for linearly coupled damped 
oscillators with damping constants ji, eigenfrequencies Wj, and external 
forcing hi(t). The other forcing terms on the right-hand side are due to 
interactions of suppliers. They appear only if J is not of diagonal, but of 
Jordan normal form with some Jis+i ^ 0. Because of hi = J^i+i, Eqs. 
(|45[l can again be analytically solved in a recursive way, starting with the 
highest index i = u. Note that, in the case D = E (i.e. dij = Sij), J a are 
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Supply Networks: From the Bullwhip Effect to Business Cycles 13 

the eigenvalues of the input matrix C and < \J%%\ < 1. Equation l|45|l has 
a special periodic solution of the form 

Xi {t) =x1e l{at -^\ (47) 
h i (t) = h° i e iat , (48) 

where i = y— T denotes the imaginary unit. Inserting this into l|45|l and 
dividing by e lat immediately gives 

{-a 2 + 2ia 7i + Lu^x ^-^ = b % {V + iaW>? +1 e- ix<+1 + fo? . (49) 

With c ±1<;4 = cos((f>) ± isin(</>) this implies 



x°e" iXl 



where 



bi(V + iaW)x° i+1 e- i ^+ 1 + ti\ ^Re 2 + Im 2 e'« 



-a 2 + 2ia 7i + u t 2 ^(w 4 2 - a 2 ) 2 + (2a 7l ) 2 e^ ' 

(50) 

Re = biX° +1 [V cos(xi+i) + aW sin(xi+i)] + 

Im = biX° +1 [aW cos(xi+i) - V sm(xi+i)} , (51) 



and 



with 



l [V 2 + (aW) 2 ](b l x'i +1 ) 2 + h<iH l + (h^ 
X *~V (^ 2 -a 2 ) 2 + (2a 74 ) 2 



Hi = 2b iX ° i+1 [V cos(xi+i) + aWsin(xi+i)] • (53) 
Finally, we have 

Xi = <Pt-pi with tan^i = 2a 7i /(w. i 2 - a 2 ) (54) 

and 

= kx Q l+1 [aW cos{xr+i) - V"sin(xz+i)] 
1 M° +1 [Vcos( Xi+1 ) + aWsin( Xl+1 )]+/ l o- W 

For ft,? = 0, we obtain 

tan((^ 4 - Xi) = tan((5 - x-i+i) (56) 

with 

tan 5 = aW/V , (57) 
i.e. the phase shift between i and i + 1 is just 

Xi ~ Xi+i =<Pi-5. (58) 



February 2, 2008 11:33 



WSPC/Trim Size: 9in x 6in for Review Volume 



hclbinglaemmer 



14 Dirk Helbing and Stefan hammer 

According to Eq. i|45f) . the dynamics of our supply network model can be 
reduced to the dynamics of a sequential supply chain. However, the eigen- 
values are now potentially complex and the new enti ties i h ave the meaning 
of ' 'quasi- suppliers" (analogously to "quasi-species'®E3) defined by the 
linear combination x(t) = T _1 <f(r). This transformation makes it possi- 
ble to define the bullwhip effect for arbitrary supply networks: It occurs 
if the amplitude x® is greater than the oscillation amplitude x® +1 of the 
next downstream supplier i + 1 and greater than the amplitude h® of the 
external forcing, i.e. if x°/max(x° +1 , h°) > 1. 

Note that in contrast to sequential supply chains, the oscillation am- 
plitude of x(t) may be amplified in the course of time, depending on the 
network structure. This case of absolute instability can occur if at least 
one of the eigenvalues A^± of the homogeneous equation l|45ll resulting for 
hi = bi = has a positive real part, which may be true when some complex 
eigenvalues J a exist. The (up to) 2u eigenvalues 

= -[l + W(l- J„)]/2 ± y/[l + W(l - J l4 )] 2 /4 - V(l - .hi) (59) 

depend on the (quasi-)supplier i and determine the temporal evolution of 
the amplitude of deviations from the stationary solution. 

3.2. Examples of Supply Networks 

It is useful to distinguish the following cases: 

Symmetric supply networks: If the supply matrix S is symmetric, as 
for most mechanical or electrical oscillator networks, all eigenvalues J a are 
real. Consequently, if u>i < ji (i.e. if V is small enough), the eigenvalues 
\i t ± of M are real and negative, corresponding to an overdamped behavior. 
However, if W is too small, the system behavior may be characterized by 
damped oscillations. 

Irregular supply networks: Most natural and man-made supply net- 
works have directed links, and S is not symmetric. Therefore, some of the 
eigenvalues Ja will normally be complex, and an overdamped behavior is 
untypical. The characteristic behavior is rather of oscillatory nature (al- 
though asymmetry does not always imply complex eigenvalues^, see Fig.^J. 
For small values of W, it can even happen that the real part of an eigen- 
value Aj,± becomes positive. This implies an amplification of oscillations in 
time (until the oscillation amplitude is limited by non-linear terms). Sur- 
prisingly, this also applies to most upper triangular matrices, i.e. when no 
loops in the material flows exist. 
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Regular supply networks: Another relevant case are regular supply 
networks. These are mostly characterized by degenerate zero eigenvalues 
Ju = and Jordan normal forms J, i.e. the existence of non-vanishing 
upper-diagonal elements Ji^+\. Not only sequential supply chains, but also 
fully connected graphs, regular supply ladders, and regular distribution sys- 
tems belong to this case^l (see Fig. 2^, c > d). This is characterized by the 
two M-fold degenerate eigenvalues 



A± = -(1 + W)/2 ± y/(l + W) 2 /A - V , (60) 

independently of the suppliers i. For small enough values V < (1 + W) 2 /4, 
the corresponding supply systems show overdamped behavior, otherwise 
damped oscillations. 



Fig. 4. Different examples of supply networks: (a) Sequential supply chain, (b) closed 
supply chain or supply circle, (c) regular supply ladder, (d) regular hierarchical distri- 
bution network. 



For the purpose of illustration, the following equations display some 
regular input matrices C and their corresponding Jordan matrices J: For a 
fully connected network we have 
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The Jordan normal matrix J of a sequential supply chain corresponds to 
the input matrix C itself, i.e. J = C. This, however, is quite exceptional. 
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For the supply ladder displayed in Fig. we have 
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(62) 



/0 10000000 0\ 
0010000000 
0001000000 
0000100000 
0000000000 
0000000000 
0000000000 
0000000000 
0000000000 

Vo 00000000 o/ 



where the number of ones corresponds to the number of levels of the supply 
ladder. For the hierarchical distribution network shown in Fig.^Ji, but with 
3 levels only, we have 



C = 



1 i 

2 2 



0\ 

1 1 



ii 

\ \ 




Vo 000000/ 



and J 



fo 1 0\ 
10 

10 

1 

\0 000000/ 



(63) 



Note that the Jordan normal forms of different input matrices C may be 
identical, but the transformation matrix T and the driving term h(t) would 
then be different. 

Randomized regular supply networks belong to the class of irregu- 
lar supply networks, but they can be viewed as slightly perturbed regular 
supply networks. For this reason, there exist approximate analytical re- 
sults for their eigenvalues. Even very small perburbations of the regular 
matrices S discussed in the previous paragraph can change the eigenvalue 
spectrum qualitatively. Instead of the two multiply degenerate eigenvalues 
A± of Eq. f^|l , we find a scattering of eigenvalues around these values. The 
question is: why? 

In order to assess the behavior of randomized regular su pply networks, 
we apply Gersgorin's theorem on the location of eigenvalues!^ According 
to this, the n complex eigenvalues Afc S C of some n x n-matrix N are 
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located in the union of n disks: 

A fe g [J \ z G C : \z - N u \ < ^ \ (64) 

Furthermore, if a union of I of these n discs form a connected region that 
is disjoint from all the remaining n — I discs, then tere are precisely / eigen- 
values of N in this region. 

Let us now apply this theorem to the perturbed matrix 

C„ = C + 77P . (65) 

For small enough values of 77, the corresponding eigenvalues Ja should be 
located within discs of radius 

R i (v)=vJ2\ p i? ) \ ( 66 ) 

3 

around the (possibly degenerated) eigenvalues Ja 0/ the original matrix C 
(see Fig- - This radius grows monotonously, but not necessarily linearly 
in the parameter 77 with < i] < 1, which allows to control the size of the 
perturbation. Moreover, = R^PR,,, where R ); is the orthogonal ma- 
trix which transforms C v to a diagonal matrix D™ , i.e. R^ 1 C I) R^ = D«' . 
(This assumes a perturbed matrix C v with no degenerate eigenvalues.) Sim- 
ilar discs as for the eigenvalues of C v can be determined for the associated 
eigenvalues of the perturbed (2u x 2u)-matrix M, belonging to the 
perturbed u x u-matrix see Eq. (|38|l and Fig. [S] 



Fig. 5. (a) Example of a supply network (regular distribution network) with two degen- 
erated real eigenvalues (see squares in the right subfigure). (b) The eigenvalues of 
the randomly perturbed supply network (crosses) are mostly complex and located within 
Gersgorin's discs (large circles). (After Refill) 

Let us now discuss the example of a structural perturbation of a se- 
quential supply chain (with 7/ = 0) towards a supply circle (with 77 = 1), 
see Fig. 0Jd. For this, we set 
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(67) 
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While the normal form for r\ — is given by a Jordan matrix J which agrees 
with C, for any rj > we find the diagonal matrix 



/Jn ... \ 
J 22 • • • 



(68) 



V J nn ) 

where the diagonal elements Jn are complex and equally distributed on a 
circle of radius ^/rj around the origin of the complex plane. Therefore, even 
an arbitrarily small perturbation can change the eigenvalues qualitatively 
and remove the degeneration of the eigenvalues. 

4. Network-Induced Business Cycles (with Ulrich Witt and 
Thomas Brenner) 

In order to investigate the macroeconomic dynamics of national economies, 
it is useful to identify the production units with economic sectors (see 
Fig. El . It is also necessary to extend the previous supply network model 
by price-related equations. For example, increased prices Pi(t) of products 
of sector i have a negative impact on the consumption rate Yi (t) and vice 
versa. We will describe this by a standard demand function Li with a neg- 
ative derivative L'^Pi) = dLi(Pi) / dPf. 

Y i (t) = [Y° + Z i (t)]L i (P i (t)). (69) 

This formula takes into account random fluctuations £j (t) over time around 
a certain average consumption rate 5^° and assumes that the average value 
of Lj(Pj(i)) is normalized to one. The fluctuation term £j(t) is introduced 
here in order to indicate that the variation of the consumption rate is a 
potentially relevant source of fluctuations. 



Fig. 6. Main service and commodity flows among different economic sectors according 
to averaged input-output data of France, Germany, Japan, UK, and USA. For clarity of 
the network structure, we have omitted the sector 'wholesale and retail trade', which is 
basically connected with all other sectors. 



Inserting 169fl into results in 
dNi(t) 



Qi(t) - E ^C*) - K° + &(t)]£i(Pi(*)) . (70) 



dt 

J_1 =n(t) 
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Herei n, w e have assumed = <5y as in Leontief's classical input-output 
model!™ Moreover, in our simulations we have applied the common linear 
demand function 

Li(Pi) = m S x{0,L° i -L 1 i P i ), (71) 

where and L\ are non-negative parameters. 

Due to the price dependence of consumption, economic systems have the 
important equilibration mechanism of price adjustment. These can compen- 
sate for undesired inventory levels and inventory changes. A mathematical 
equation reflecting this is 

1 dPi ( N? \ m dNj (?2) 



Pi{t) dt \Ni(t) J Ni(t) dt 

The use of relative changes guarantees the required non-negativity of prices 
Pi{t) > 0. Vi is an adaptation rate describing the sensitivity to relative 
deviations of the actual inventory iVj(t) from the desired one Nf, and fit is a 
dimensionless parameter reflecting the responsiveness to relative deviations 
(dNi j dt) I Ni(t) from the stationary equilibrium state. 

If the same criteria arc applied to adjustments of the production rates 
Qi(t), we have the equation 

1 dQ % _ . ( Nf \ _ Ji^m 



Qi(t) dt \Ni(t) J Ni(t) dt 

oti = Vi/vi is the ratio between the adjustment rate of the output flow and 
the adjustment rate of the price in sector i. For simplicity, the same ratio 
oti — fri/ Hi will be assumed for the responsiveness. 

4.1. Treating Producers Analogous to Consumers (with 
Dieter Armbruster) 

According to Eqs. I|tj9|) an d l|71|l . the change of the consumption rate in 
time is basically given by 

f-W+M^f-W+iWW?, c«> 

i.e. it is basically proportional to the price change dPi/dt, but with a neg- 
ative and potentially fluctuating prefactor. However, 

1 dQi 1 dP t 

TTTTT - TT ' ( 75 ) 



Qi(t) dt Pi(t) dt 

i.e. according to Eqs. 172JI and J22J, the change dQi/dt of the production 
rate close to the equilibrium state with Qi(t) w Q° and Pi(t) w Pf is 
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basically proportional to the change dPi/dt in the price with a positive 
prefactor. Is this an inconsistency of the model? Shouldn't we better treat 
producers analogous to consumers? 

In order to do so, we have to introduce the delivery flows D^j of products 
i to producers j. As in Eq. l(7i|. we will assume that their change dDij/dt 
in time is proportional to the change dPi/dt in the price, with a negative 
(and potentially fluctuating) prefactor. We have now to introduce the stock 
level Oi (t) in the output buffer of producer i and the stock levels (t) of 
product i in the input buffers of producer j. For the output buffer, we find 
the balance equation 

^ = Qi{t)- y £D ij (t)-Y i (t) (76) 

3 

analogous to Eq. 1)70(1 . as the buffer is filled with the production rate of 
producer i, but is emptied by the consumption rate and delivery flows. For 
the input buffers we have the balance equations 

^-=D ij {t)-c ij Q j {t), (77) 

as these are filled by the delivery flows, but emptied with a rate proportional 
to the production rate of producer j, where c^- are again the input coeffi- 
cients specifying the relative quantities of required inputs i for production. 
Generalizations of these equations are discussed elsewherePSI 

Let us now investigate the differential equation for the change dNi/dt 
of the overall stock level of product i in all input and output buffers. We 
easily obtain 

= + ? dJ i = Qi ® E - m , m 

3 3 

as before. That is, the equations for the delivery flows drop out, and we stay 
with the previous set of equations for Ni, Qi, Pi, and lj. As a consequence, 
we will focus on Eqs. I|69|) to (|73|l in the following. 



5. Reproduction of Some Empirically Observed Features of 
Business Cycles 

Our simulations of the above dis-aggregate (i.e. sector-wise) model of 
macroeconomic dynamics typically shows asynchronous oscillations, which 
seems to be characteristic for economic systems. Due to phase shifts be- 
tween sectoral oscillations, the aggregate behavior displays slow variations 
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of small amplitude (see Fig.0). If the function Lj(Pj) and the parameters 
Vil Hi 2 are suitably specified, the non-linearities in Eqs. 170(1 to l(69[) will 
additionally limit the oscillation amplitudes. 



Fig. 7. Typical simulation result of the time-dependent gross domestic product 
EjQiM-PiM m percent, i.e. relative to the initial value. The input matrix was cho- 
sen as in Figs.E]and|nii-d, but Y? was determined from averaged input-output data. Q' 
was obtained from the equilibrium condition, and the fluctuations £;(t) were specified 
as a Gaussian white noise with mean value and variance a = 10.000 (about 10% of 
the average final consumption). The initial prices P;(0) were selected from the interval 
[0.9;1.1]. Moreover, in this example we have assumed Li(Pi) = max[0, 1 + d(Pi — P-)] 
with d = f'(P?) = -L\ = -10 and the parameters v t = 0.1, m = 0.0001, a l = 1 = P9, 
and TV? = . Although this implies a growth of small oscillations (cf. Fig. Id), the 
oscillation amplitudes are rather limited. This is due to the non-linearity of model equa- 
tions 1701 to 1691 and due to the phase shifts between oscillations of different economic 
sectors i. Note that irregular oscillations with frequencies between 4 to 6 years and am- 
plitudes of about 2.5% are qualitatively well compatible with empirical business cycles. 
Our material flow model can explain uj-shaped, non-periodic oscillations without hav- 
ing to assume technological shocks or externally induced perturbations. The long-term 
growth of national economies was intentionally not included in the model in order to 
separate this effect from network-induced instability effects. (After RefEJ) 

Our business cycle theory differs from the dominating one^E^ in sev- 
eral favourable aspects: 

(i) Our theory explains irregular, i.e. non-periodic oscillations in a 
natural way (see Fig. EJ. For example, iw-shaped oscillations result 
as superposition of the asynchronous oscillations in the different 
economic sectors, while other theories have to explain this observa- 
tion by assuming external perturbations (e.g. due to technological 
innovations). 

(ii) Although our model may be extended by variables such as the labor 
market, interest rates, etc., we consider it as a potential advantage 
that we did not have to couple variables in our model which are 
qualitatively that different. Our model rather focusses on the ma- 
terial flows among different sectors. 

(iii) Moreover, we will see that our model can explain emergent oscilla- 
tions, which are not triggered by external shocks. 



February 2, 2008 11:33 



WSPC/Trim Size: 9in x 6in for Review Volume 



hclbinglaemmer 



22 



Dirk Helbing and Stefan hammer 



Fig. 8. The displayed phase diagram shows which dynamic behavior is expected by our 
dis-aggregate model of macroeconomic dynamics depending on the respective parameter 
combinations. The individual curves for different countries are a result of the different 
structures of their respective input matrices C. As a consequence, structural policies can 
influence the stability and dynamics of economic systems. 

5.1. Dynamic Behaviors and Stability Thresholds 

The possible dynamic behaviors of the resulting dis-aggregate macroeo- 
nomic model can be studied by analytical investigation of limiting cases 
and by means of a linear stability analysis around the equilibrium state 



(in which we have N^t) = Nf, Y t (t) = Y?, Q° - EjCyQj = Y?, and 
P%(t) = Pj°). The linearized equations for the deviations Jij(t) = Ni(t) — N® , 
Pi(t) = Pi(t) — Ff, and qi{t) — Qi(t) — Q° from the equilibrium state read 



This system of coupled differential equations describes the response of the 
inventories, prices, and production rates to variations in the demand. 
Denoting the m eigenvalues of the input matrix C = (cy) by J a with 
| Jit | < 1, the 3m eigenvalues of the linearized model equations are (m 
times) and 



dni 
~dt 



(79) 



3 




(80) 



(81) 




(82) 



where 



Ai = m[d + diA(l - Ju)} , 
Bi = Vi[Ci + &iA(l - J«)] ) 

Q = ppy i °|/?(pP)l/^ ; 

A = Q°i/N? , 



(83) 



see Fig. |H1 Formula 182(1 becomes exact when the matrix C is diagonal 
or the parameters /XjCj, aj/XjA) v i@% an d cti^iDi are sector-independent 
constants, otherwise the eigenvalues must be numerically determined. 
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Fig. 9. Properties of our dynamic model of supply networks for a characteristic input 
matrix specified as average input matrix of macroeconomic commodity flows of several 
countries (top) and for a synthetic input matrix generated by random changes of input 
matrix entries until the number of complex eigenvalues was eventually reduced to zero 
(bottom). Subfigures (a), (e) illustrate the color-coded input matrices A, (b), (f) the 
corresponding network structures, when only the strongest links (commodity flows) are 
shown, (c), (g) the eigenvalues Ja = He(Ja) + ilm(Ja) of the respective input matrix A, 
and (d), (h) the phase diagrams indicating the stability behavior of the model equations 
1701 to 1691 on a double-logarithmic scale as a function of the model parameters ct; = a 
and Ui/fii 2 = u/fj, 2 = V/M 2 . The other model parameters were set to in = C; = Di = 
P9 = iV? = YP = 1. Surprisingly, for empirical input matrices A, one never finds an 
overdamped, exponential relaxation to the stationary equilibrium state, but network- 
induced oscillations due to complex eigenvalues J;;. (After Ref@) 

It turns out that the dynamic behavior mainly depends on the parame- 
ters a,-, 2/j/ Hi 2 , and the eigenvalues J„ of the input matrix C (see Fig. 
In the case &i — > of fast price adjustment, the eigenvalues \ t ± are given 



by 



2\ h± = -md ± y/imCi)* - 4viQ , (84) 



i.e. the network structure does not matter at all. We expect an exponential 
relaxation to the stationary equilibrium for < Vij ft 2 < Ci/A, other- 
wise damped oscillations. Therefore, immediate price adjustments or simi- 
lar mechanisms are an efficient way to stabilize economic and other supply 
systems. However, any delay (dj > 0) will cause damped or growing oscil- 
lations, if complex eigenvalues Ja — Ke(Ju) +iIm(J„) exist. Note that this 
is the normal case, as typical supply networks in natural and man-made 
systems are characterised by complex eigenvalues (see top of Fig. |5J. 
Damped oscillations can be shown to result if all values 

i>i/{J,i 2 = a^i/fi 2 (85) 

lie below the instability lines 

Villi 2 » \Ci + &iDi[l - Re(J«)]} 

v 1 + pMfef ) (86) 

given by the condition Re(Aj i ±) < 0. For identical parameters vij [i 2 — 
v/fj 2 and on — a, the minimum of these lines agrees exactly with the 
numerically obtained curve in Fig. [5] Values above this line cause small 
oscillations to grow over time (see Appendix A) . 
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In some cases, all eigenvalues Ju of the input matrix C are real. This 
again applies to symmetric matrices C and matrices equivalent to Jordan 
normal forms. Hence, the existence of loops in supply networks is no suffi- 
cient condition for complex eigenvalues Ju (see also Fig. If). It is also no 
necessary condition. For cases with real eigenvalues only, Eq. 1)82(1 predicts 
a stable, overdamped behavior if all values Vi/fa — a^i/ fii 2 lie below the 
lines 

Vi/Hi 2 « [d + OiA(l - J«)]/4 (87) 

defined by min^A^ 2 — 4Bj) > (see Appendix B). For identical parameters 
Vi//Ji 2 = v I /i 2 and en — a, the minimum of these lines corresponds ex- 
actly to the numerically determined curve in Fig. El Above it, one observes 
damped oscillations around the equilibrium state, but growing oscillations 
are not possible. In supply systems with a slow or non-existent price ad- 
justment mechanism (i.e. for on 1 or d = 0), Eq. 1(87(1 predicts an 
overdamped behavior for real eigenvalues Ju and 

Ui/fii 2 < A(l - Ju)/4, (88) 

for all i, while Eq. ((86(1 implies the stability condition 

hi fa 2 < A[l - He(J«)]{l + [1 - Re(J u )} 2 /lm(J u ) 2 } (89) 

for all i, given that some eigenvalues Ju are complex. Moreover, for the 
case of sector-independent constants V = a^Di and W = &iPiDi, the 
eigenvalues A^± can be calculated as 

K± = -W(l - ± y/[W(l - J i4 )] 2 /4 - V(l - J u ) ■ (90) 

Considering that Eq. I(3(j|) for the adjustment of production speeds is slightly 
more general than specification 1(81(1 . this result is fully consistent with 

Eq. ijsnj)- 

6. Summary 

In this contribution, we have investigated a model of supply networks. It is 
based on conservation equations describing the storage and flow of inven- 
tories by a dynamic variant of Leontief 's classical input-output model. A 
second set of equations reflects the delayed adaptation of the production 
rate to some inventory-dependent desired production rate. Increasing values 
of the adaptation time T tend to destabilize the system, while increasing 
values of the time horizon At over which inventories are forecasted have a 
stabilizing effect. 
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A linear stability analysis shows that supply networks are expected to 
show an overdamped or damped oscillatory behavior, when the supply ma- 
trix is symmetrical. Some regular supply networks such as sequential supply 
chains or regular supply hierarchies with identical parameters show stable 
system behavior for all values of the adaptation time T. Nevertheless, one 
can find the bullwhip effect under certain circumstances, i.e. the oscillations 
in the production rates are larger than the variations in the consumption 
rate (demand). The underlying mechanism is a resonance effect. 

Regular supply networks are characterized by multiply degenerate eigen- 
values. When the related supply matrices are slightly perturbed, the de- 
generation is broken up. Instead of two degenerate eigenvalues, one finds 
complex eigenvalues which approximately lie on a circle around them in 
the complex plane. The rela ted heterogeneity in the model parameters can 
have a stabilizing effect P^l when they manage to reduce the resonance ef- 
fect. Moreover, a randomly perturbed regular supply network may become 
linearly unstable, if the perturbation is large. If the supply matrix is irreg- 
ular, the supply network is generally expected to show either damped or 
growing oscillations. In any case, supply networks can be mapped onto a 
generalized sequential supply chain, which allows one to define the bullwhip 
effect for supply networks. 

While previous studies have focu ssed o n the synchronization of oscil- 
lators in different network topologies j ^ l ^H we have found that many sup- 
ply networks display damped oscillations, even when their units — and linear 
chains of these units — behave in an overdamped way. Furthermore, networks 
of damped oscillators tend to produce growing ( and mostly asynchronous ) 
oscillations. Based on these findings, it is promsing to explain business cy- 
cles as a result of material flows among economic sectors, which adjust their 
production rates in a decentralized way. Such a model can be also extended 
by aspects like labor force, money flows, information flows, etc. 

7. Future Research Directions 
7.1. Network Engineering 

Due to the sensitivity of supply systems to their network structure, network 
theory^ ^ l 1 M-*- H 1 ^ 1 1 3 1 5^ 1 can raake useful contributions: On the basis of 
Eqs. I|86(l and 187(1 one can design stable, robust, and adaptive supply net- 
works ( "network engineering" ) . For this reason, our present studies explore 
the characteristic features of certain types of networks: random networks, 
hierarchical networks, small world networks, preferential attachment net- 
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works, and others. In this way, we want to identify structural measures 
which can increase the robustness and reliability of supply networks, also 
their failure and attack tolerance. These results should, for example, be 
relevant for the optimization of disaster management 

7.2. Cyclic Dynamics in Biological Systems 

Oscillations are probably not always a bad feature of supply systems. For 
example, in s yste ms with competing goals (such as pedestrian counterflows 
at bottleneck]^] or intersecting traffic streams), oscillatory solutions can 
be favourable. Oscil lations are a lso found in ecological systemJ^ and some 
metabolic processes ESEZESE2l w hi c h could be also treated by a generalized 
model of supply networks. Examples are oscillations of nutrient f low s in 
slime molds the citrate cy cle !^ the glycolytic oscillations in yeast and 
the Calcium cycle?S2S2ISI Considering the millions of years of evolution, 
it is highly unlikely that these oscillations do not serve certain biological 
functions. Apar t fro m metabolic flows, however, we should also mention 
protein networkJ^I and metabolic expression networkJ^EHl as interesting 
areas of application of (generalized) supply network models. 

7.3. Heterogeneity in Production Networks 

Another interesting subject of investigation wou ld b e heterogeneity. The 
relevance of this issue is known from traffic flowd^ and multi-agent sys- 
tems in economical Preliminary results indicate that, depending on the 
network structure, heterogeneous paramet ers of production units can re- 
duce the bullwhip effect in supply networks!^ That is, identical production 
processes all over the world may destabilize economic systems (as mono- 
cultures do in agriculture). However, in heterogeneous systems a lot more 
frequencies are influencing a single production step, which may also have 
undesireable effects. These issues deserve closer investigation. 

7.4. Multi-Goal Control 

Let us go back to Eq. 1)34(1 for the control of supply or production rates. 
This equation is to be applied to cases where a production unit produces 
several products and the production rate is not only adapted to the inven- 
tory of one dominating product, but also to the inventories of some other 
products. This leads to the problem of multi-goal control. It is obvious that 
the different goals do not always have to be in harmony with each other: 
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A low inventory of one product may call for an increase in the production 
rate, but the warehouse for another product, which is generated at the same 
time, may be already full. A good example is the production of chemicals. 
We conjecture that multi-goal control of production networks may imply 
additional sources of dynamical instabilities and that some of the problems 
may be solved by price adjustments, which can influence the consumption 
rates. Again, a closer investigation of these questions is needed. 

7.5. Non- Linear Dynamics and Scarcity of Resources 

Finally, we should note that coupled equations for damped harmonic oscil- 
lators approximate the dynamic behavior of supply networks only close to 
their stationary state. When the oscillation amplitudes become too large, 
non-linearities will start to dominate the dynamic behavior. For example, 
they may select w ave mo des or influence their phases in favour of a syn- 
chronized behavio r P" 1 _ ' etc . Deterministically chaotic behavior seems to 
be possible as well P^^^^^ Numerical results for limited buffer sizes 
a nd t ransport capacities have, for example, been presented by Peters et 



Additional non-linearities come into play, if there is a scarcity of re- 
sources required to complete certain products. Let's assume that the rate 
of transfering products i to j is Zj+ . If iVj is the number of products that may 
be delivered, the maximum delivery rate is therefore limited by ZjiNi(t). 
The maximum production rate is given by the minimum of the delivery rate 
of all components i, divided by the respective number cy of parts required 
to complete one unit of product i. In mathematical terms, we have to make 
the replacement 



in Eq. (|70|l and corresponding replacements in all derived formulas!*^ 
Therefore, the generalized relationship for systems in which resources may 
run short is 



a m 




(91) 



dt 



^2{dij - Cij)Qj{t) mm ( 1, 



Z ]k N k {t) 
c jk 



) 



(92) 



These coupled non-linear equations are expected to result in a complex 
dynamics, if the transportation rate Zj k or inventory N k are too small. 
Note that this non-linearity may be even relevant close to the stationary 
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state. In such cases, results of linear stability analyses for systems with 
scarce resources would be potentially misleading. 

In situations where scarce resources may be partially substituted by 
other available resources, one may use the replacement 



Qj(t) 



. 1 1 / z i 



(93) 



instead of l|91l) . Such an approach is, for example, reasonable for disaster 
management P22z.,- is a parameter allowing to fit the ease of substitution of 
resources. The minimum function results in the limit Zj — > — oo. 
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Appendix A. Boundary between Damped and Growing 
Oscillations 

Starting with Eq. Ij82(l . stability requires the real parts Re(Aj) of all eigen- 
values Xi to be non-positive. Therefore, the stability boundary is given by 
maxi Re(Aj) =0. Writing 

C i + a i D i (l-J ii ) = e i + i (A.l) 

with C t = P?Y?\fl(P?)\/N? and defining 

§i = C i +a i D i [l-Re(J u }} ) 

pi = ^faiDilm(Ja) (complex conjugate eigenvalues), 

7i = , (A.2) 

we find 



2A»//i 2 = -§i - ifa + \/R t + Hi (A.3) 

with 

84=8? and I i =2Q i h-lih- (A.4) 
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The real part of l|A.3|) can be calculated via the relation 



Bc(y/Ri ± i/i) = \/1(VR? + h 2 + Rt) • (A.5) 
The condition Re(2Ai//x,) = is fulfilled by % — and 

= 4Si(l + 6i 2 /fo 2 ) , (A.6) 
i.e. the stable regime is given by 

7i Vi 



2 

j 



for all i, corresponding to Eq. (|86|l .^ 



Appendix B. Boundary between Damped Oscillations and 
Overdamped Behavior 

For &i > 0, the imaginary parts of all eigenvalues vanish if lm(Ju) — 
(i.e. 0i = 0) and if Ri > 0. This requires 

^ = li < k - = §i = Ci + &iA(l - J«) (B.l) 

for all i, corresponding to Eq. I|87|l .^ 
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